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I.  INTRODUCTION 

A most  common  technique  in  data  analysis  involves  fitting  a function, 
the  form  of  which  is  dictated  by  theoreticiil  considerations,  to  experimental 
data.  The  parameters  of  the  function,  which  are  adjusted  to  provide  the 
best  fit  to  the  data,  are  related  to  the  physical  quantities  which  the 
experiment  sought  to  measure;  hence,  the  values  of  these  parameters 
become  the  desired  results  of  the  experiment.  In  particular,  most  of  the 
experiments  conducted  at  this  Laboratory,  including  recent  work  on  energy 
and  rate  dependence  of  solid  state  x-ray  dosimeters y response  of 
ferroelectric  dosimeters,  yield  outputs  which  vary  with  changes  in 
experimental  conditions.  By  a theoretical  analysis  of  the  mechanisms 
involved  in  the  interaction  of  radiation  with  these  dosimeters,  their 
expected  responses  have  been  described  as  functions  of  the  experimental 
conditions.  From  the  parameters  of  these  functions,  the  energy,  dose 
and  rate  dependence,  time  dependence,  electrical  characteristics, 
have  been  determined.  In  some  cases  the  failure  of  the  expected 
function  to  satisfactorily  fit  the  data  has  resulted  in  the  discovery 
of  other  mechanisms  involved  in  the  response,  and  in  the  subsequent 
generation  of  a more  complete  functional  description. 

Because  of  the  broad  application  of  this  kind  of  analysis,  data 
fitting  routines  are  commonplace;  most,  however,  involve  fitting  a 
particular  functional  form,  such  as  a power  series  or  Fourier  series. 

For  experiments  whose  functional  behaviour  is  not  such  a series,  the 
subsequent  extraction  of  t-he  desired  quantities  may  be  very  difficult. 

FNFIT  is  designed  to  allow  the  analyst  to  input  his  own  functional 
form,  and  reap  those  parameters  which  are  directly  applicable  to  the 
experiment.  The  program  was  intended  to  be  quick  and  easy-to-use  for 
the  analyst  who  has  only  occasional  need  to  fit  a simple  function, 
yet  versatile  enough  ’.o  handle  coi.plicated  fu^-rtions,  and  functions  which 
change  erratically  wi :h  small  ''hanges  in  >arameters.  The  user  supplied 
function,  G (x;  A(l),  A(2)...)  can  contain  up  to  20  parameters,  A(J) , and 
up  to  1000  data  points,  [X(I),  Y(lj],  can  be  input. 
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Five  or  six  sets  (the  sixth  being  optional)  of  cards  are  read  in  to 
run  the  program.  These  sets  input  the  following: 

EXEC  A.  The  function,  G. 

EXEC  B.  Assemble  and  execute  instructions 
SET  A.  Control  information 
SET  B.  Starting  values  for  the  parameters 
SET  C.  Limits  on  the  parameters  (Optional) 

SET  D.  The  Data 

These  card  sets  are  described,  along  with  some  "convenience"  options, 
are  noted  in  the  following  section,  ("PROGRAM  INPUT"). 

The  program,  at  the  end  of  the  run,  outputs  the  parameters,  the 
data  points  and  function  value  at  each  point,  and  the  derivative  of  the 
function  at  each  point.  Various  outputs  are  available  at  the  beginning 
of  each  iteration. 


II.  PROGRAM  INPUT 
A.  Blank  (Frce-Field)  Format 

Since  FNFIT  makes  use  of  the  BLANK  FORMAT  feature  of  the  1108, 
a note  of  description  is  in  order.  In  BLANK  format,  the  computer, 
directed  by  a READ  statement,  begins  to  scan  a record  (card), 
interpreting  whatever  it  finds  as  the  type  of  variable  sought  by  the 
statement;  the  field  ends  at  a comma  or  at  the  end  of  the  record.  Blanks 
are  usually  ignored,  except  an  entirely  blank  field  is  read  as  integer 
zero. 

For  example:  READ  (5,9)  A,  B,  J 

9 FORMAT  ( ) 

would  read 

9.0  , 6.7E-1,7 

or 

9.0  , 6.7-1 

7 


6 


5S.y^  ifK 


I 

F 

t 

i 

I 

; 


identically,  giving  the  valuer. 

A = 9.0,  B = .67,  and  J - 7 
The  advantages  of  using  BLANK  Format  are: 

1.  Numbers  need  not  be  carefully  placed  in  fields 

on  the  cards.  One  merely  punches  the  numbers,  separated  by  commas. 

2.  The  same  read  and  format  statement  can  read  data  on  one  or 
several  cards,  decided  by  the  user  at  data  input  time.  This  is  useful 
in  FNFIT. 

One  precaution  to  remember  is  that  the  end  of  a card  acts  like  a comma 
(unless  a comma  is  in  column  80).  Hence  the  data 

9.0  , 6.7E-1, 

7 

would  set  A = 9.0,  B = .67,  J = 0,  and 
read  7 on  next  READ  statement. 

B.  Program  Operation 

The  six  card  sets,  described  in  the  introduction,  are  constructed  as 
follows.  (Control  cards  (beginning  with  0)  are  for  the  EXEC  VIII 
operating  system  for  the  UNIVAC  1108). 

EXEC  A 

The  first  set  of  cards  must  cause  Subroutine  FCN  to  be  compiled 
with  the  function  G inserted  where  indicated  [See  Listing  - Appendix  A]. 

The  exact  commands  depend  on  the  file  names  used,  etc.  At  the  APG-EA 
UNIVAC  1108  facility,  this  requires  only  an  executive  card, 

0ADD  AMUCK*FNFIT.D0-G 

followed  by  FORTRAN  statements  that  construct  the  user  supplied  function, 

G(x;  A(l),  A(2),  ,..).  For  example 

G = A(l)  + A(2)  * X **  2 + X **  A(3) 
or 

G = 0. 
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DO  10  I = 1,  4 


10  G = A(I)  * SIN  (I  * A(5)  ♦ X)  + G 

The  0ADD  conunand  enters  the  instructions  from  element  DO-G  ir.co  the 


run  stream. 


EXEC  B 


The  next  instructions  must  make  aji  absolute  element  which  includes 
the  newly  compiled  and  execute  the  program.  Again^  exact  commands 
depend  on  the  system  which  in  our  case  requires  only  one  card  - 
0ADD  AMUCK*FNFIT.DO-FNFIT 

The  program  is  now  running.  The  following  sets  insert  data  for  a 
particular  fit  using  the  user-supplied-function. 

SET  A 

(1)  ‘’Gab  cards"  - i.e.  Any  cards  which  do  not  start  with  five  integers 
separated  by  commas,  will  be  printed  at  the  top  of  any  output.  This  allows 
printing  the  G-function  or  messages. 

The  first  cards  starting  with  five  integers  will  be  read  as  the 
following  control  card. 

(2)  "Control  card"  ( FORMAT  ( ) ) 


MXSTP,  MAXIT,  IPRT,  ICHK,  NHOLD 


where 


MXSTP:  Maximum  number  of  steepest  gradient  steps 

MAXIT:  Maximum  number  of  least  square  iterations. 

IPRT:  OUTPUT  control 

= 0;  Output  only  at  end 

= 1;  Gives  parameters  and  after  every  step/iteration 
= 3;  Gives  parameters,  x^»»and  fit  to  function  after 
every  step/iteration 
= 4,  Acts  like  3,  but  includes  dG/dX 


S*  ♦ 


I 


= 2;  Acts  like  3 if  x^  has  improved,  otherwise  like  1. 
ICHK:  COOLIT  control 

= 1;  Check  parameters  for  limits  on  minima,  maxima  and 

variation  (i.e.  Call  COOLIT  each  iteration) 

= 0;  No  check 

NHOLD:  Number  of  parameters  NOT  allowed  to  vary. 

If  NHOLD  is  not  0;  then  the  next  cards  are 
(3),.,.  IHOLD(J)  ( FORMAT  ( ) ), 

the  indicies  of  the  NHOLD  parameters  which  are  to  be  held  fixed. 

For  example,  if  A(3)  is  to  be  a fixed  value,  then  NHOLD  =1,  and  the 
following  caru  is 


3 


SET  B 


(1)  NPARAM  ( FORMAT  ( ) ) - the  total  numbei  of  parameters  - A(J)  - 

in  G(x;  A(l),  ...) 

Then  NPARAM  cards,  (for  J = 1,  2,  ...  N'^ARAM) 

(2) ,...  A(J),  MAG(J)  ( FORMAT  ( ) ) 

where  A(J)  = initial  guess  for  parameter  J 

MAC(T)  = "magnitude"  of  parameter  J (i.e.  the  approximate 
power  of  10) 

e.g.  0 ^ A(l)  ^ 100.  means  MAG(l)  -•  2 

or  -0.01  ^ A(3)  ^ 0.05  means  MAG(3)  = -2 

Option  1.  If  A(JJ  - the  initial  guess  for  a parameter  - indicates  the 
order  of  magnitude  of  parameter  (J) , then  inputting  for  MAG(J)  will 
automatically  set  MAG(J)  = log  A(J).  Clearly,  to  use  this  option, 

A(J)  i 0. 


Option  2.  0EOF  in  place  of  card  (1)  on  any  search  after  the  first, 
will  compute  starting  values,  A(J) , from  previous  magnitudes. 
Eliminates  need  to  respecify  A(J)s. 


9 


t 


I 


I 


CXjUIPLE:  For  G = A(l)  * X + A(2),  input  might  be: 

2 

0.0,  1 

l.OE-2,  5 

Subsequent  fit 
0EOF 


SET  C 

If  ICHK=1  - set  C is  used. 

( COOLIT  CONTROL  SET  ) 

Subroutine  COOLIT  limits  allowed  parameter  values 

(I),...  "Coolit  Cards"  ( FOWIAT  ( ) ) 

Jl,  LI,  DUM  ( FORMAT  ( ) ) 

where 

Jl:  Subscript  of  cooled  parameter 

LI:  Kind  of  limit 
= 1 - limit  size  of  change 

= 2 > limit  minimum  value 
= 3 - limit  maximum  value 
DUM:  Value  of  the  limit 

LAST)  LAST  CARD  MUST  BE  0EOF 

Option:  Inputting  for  a coolit  card  automatically  limits  all 

parameters  to  change  by  less  than  1/2  of  their  "magnitude"  (from  SET  B) 

EXAMPLE: 

2,  1,  0.5 

6,  2,  0.0 

i 

0EOF 

This  set  of  COOLIT  Cards 

(1)  Limits  parameter  2 to  change  by  less  than  0.5  in  any  iteration. 

(2)  Limits  parameter  6 to  values  greater  than  or  equal  to  0.0. 

(3)  Limits  all  parameters  to  change  by  less  than  1/2  their  "magnitude". 
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SET  D 

1)  WT  - ( FORMAT  ( ) ) 

where  WT  is  the  power  of  Y(I)  by  which  to  weigh  each  data  pt. 

(Note:  WT  = 0.0  means  unweighed) 

2) , X(I).  Y(I)  - ( FORMAT  ()  ) 

These  are  the  data  points  to  be  fit 

Option:  If  WT  = S,  then  the  subsequent  cards  are 

2), X(I),  Y(I),  W(I)  - C FORMAT  ( ) ) 

where  W(I)  is  the  user-supplied  weight  for  each  data  point. 

(Note:  Program  normalizes  Z W(I)  = 1.0) 

Option:  If  WT  = @EOF,  the  data  and  weight  from  the  previous  run  are  used. 

Option:  If  the  first  X(I),...  = 0EOF,  the  data  from  the  previous  run 

are  used. 

LAST  CARD)  Following  the  data,  the  last  card  must  be  @EOF. 

To  end  program,  place  (another) 

@E0F 

at  end  of  deck. 


III.  HOW  FNFIT  WORKS 

For  statistical  reasons^,  it  is  assumed  that  the  function  is  a 

measure  of  the  "poorness  of  a fit",  is  given  by 

NDATA 

v2  - ’V'  W(l)  (Y(I)-F(X(I))^ 

^ ^ NDATA- (NPARAM-NHOLD) 

I = 1 


where  the  data  points,  NDATA  in  number,  are  the  set  [X(I),  Y(I)], 

F(x;  A(l),  A(2)...,  A(NPARAM))  is  the  fitting  function  containing  the 

parameters  A(J) , (NPARAM-NHOLD)  is  the  true  number  of  variable  parameters, 

and  W(I)  is  a factor  allowing  the  data  points  (X(I),  Y(I))  to  be  weighted. 

NDATA 


(FNFIT  causes 


Wfl)  = 1.) 


I = 1 
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can  uc  considered  as  a surface  in  the  NPARAM  + 1 dimensional 
space  [x^,  A(l),  A(2)...,]  The  best  fit  is  the  lowest  point  on  that 
surface  (smallest  x^) > and  the  purpose  of  FNFiT  is  to  find  the  NPARAM  + 

1 coordinates  of  that  lowest  point.  To  accomplish  this,  FNFIT  uses 
two  search  techniques 

1 . Steepest  Gradient  (£  MXSTP  steps) 

2.  x^  Minimization  MAX IT  steps) 

Consider  the  3-D  x^  suri'ace  of  Figure  1.  Suppose  one  starts  at 

point  I.  A x^  minimization  (according  to  Reference  1)  is  unstable  here 

(convex  surface).  The  best  approach  is  to  find  the  steepest  gradient 

and  step  in  that  direction.  However,  once  getting  to  IV,  the  surface 

is  concave  and  the  steepest  gradient  search  becomes  inefficient  and 

inaccurate.  Thus,  having  sensed  that  the  gradient  at  IV  was  less 

steep  than  at  III,  FNFIT  switches  to  a x^  minimization.  In  this  phase 

of  search,  the  program  approximates  x^  linearly,  and  solves  the  set  of 

NPARAM  linear  equations  to  get  the  values  of  the  parameters,  A(J) , for 

which  „ 

3A(J)  " 

FNFIT  consists  of  a main  program  (FNFIT),  10  subroutines,  and  2 
executive  instruction  elements  (DO-G  and  DO-FNFIT) . The  specific  function 
and  operation  of  the  more  involved  routines  are  discussed  below. 

FNFIT  - Main  Program 

FNFIT  serves  three  purposes.  It  acts  as  an  executive,  calling  the 
various  operations  (such  as  the  steepest  descent  search);  it  reads  all 
data  after  the  control  card;  and  it  performs  the  x^  minimization. 

Possible  confusion  about  the  program  can  be  reduced  by  distinguishing 
between  two  groupings  of  the  parameters.  The  supplied  function  (in  FCN 
via  FN)  is  a function  of  X,  containing  parameters  A(lj,  \(2) . . .A(NPARAM) . 
However,  NHOLD  of  them  may  be  excluded  by  the  user  from  being  changed 
in  value.  The  indicies  of  the  excluded  parameters  are  in  an  array 
IHOLD(l),  IH0LD(2)...IH0LD(NH0LD) . Furthermore,  as  described  in  the 
section  on  COOLIT,  some  parameters  may  be  temporarily  held  constant 
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Error  Surfacft  for  NPARAM 
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by  the  program.  There  are  NT  of  these,  and  they  go  into 
1H0LD(IH0LD(NH0LD  + 1) , . . . IHOLD(NHOLD  + NT).)  The  number  actually  being 
varied  at  any  particular  time  is  NV  = NPARAM-NHOLD-NT.  Before  each 
iteration  the  indicies  of  the  variable  parameters  are  assembled,  in 
ascending  order,  in  the  array  KL(1),  KL(2) ,. . .KL(NV) . FNFIT  uses  this 
array  to  find  the  indicies  of  the  parameters  involved  in  that  particular 
iteration  of  the  minimization  search. 

Having  read  the  input  data  and  having  returned  from  GRAD  (steepest 
gradient  search) , FNFIT  begins  the  iterative  loop.  The  steps  are  as 
follows: 

1)  It  checks  the  variable  parameters,  releasing  one  if  appropriate 
(as  described  under  COOLIT)  and  assuring  that  NV  1.  If  NV  = 0,  tha" 
iteration  is  skipped,  and  one  parameter  is  released  on  the  following 
iteration. 


2)  It  then  evaluates  the  fitting  function  - using  the  current 
values  of  the  parameters  - at  all  the  X(I)  (data  point  abscissae);  it 
calculates  x^»  outputs  as  directed  by  the  control  card. 

3)  Next,  it  sets  up  the  KL  matrix,  and  calculates  the  coefficients  for 
the  NV  simultaneous  linear  equations 

N DATA 

W(I)(Y(I)-F(X(I)))^ 

I = 1 


3AfJl 


= 0.  Recall,  y}  “ 


and  F(Z;  A(l) , A(2),...)  is  being  linearly  approximated  as 


F(Z;  A(l),  A(2),...)  = F'(Z)  + 


where  primed  quantities  refer  to  evaluation  at  tl:c  most  recent  values 
of  A(J).  Using  this,  the  NV  simultaneous  equations  can  be  written  in 
matrix  form  as 


• AA  = R 
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where 


and 


I = 1 


N DATA 

Rj  = ^ W(I)  • (Y(I)-F'(X(I)))  • 

I = 1 


/ dF(X(D) 
I dA(Kj^ 


)■ 


/ dF(X(D) 
y dA(J) 


)■ 


y 

FNFIT  gets  the  necessary  derivatives  from  DFN  (Z) . It  even  checks 
to  make  sure  that  (dF/dA(J))  is  not  identically  zero  for  some  J (i.e.,  to 
make  sure  that  F is  a function  of  all  A(J)s). 


4)  FNFIT  then  calls  LSINEQ  from  the  1108  System  Library,  and  solves 
for  the  AA  s. 

5)  At  this  point,  the  AA’s  are  checked.  If  all  are  near  zero  (no 
parameters  are  changing) , a fit  has  been  found,  and  FNFIT  exits. 

6)  Assuming  no  fit  yet,  FNFIT  goes  about  updating  the  parameters.  It 
first  calls  COOLIT  to  assu^'e  that  no  parameter  is  changed  too  much.  Upon 
return,  the  allowed  changes  are  added  to  the  appropriate  A(J)s.  COOLIT 

is  then  recalled  to  assure  that  no  parameter  has  exceeded  its  lower  or 
upper  absolute  limits. 

FNFIT  then  returns  to  step  1. 

Upon  exit,  either  via  step  5,  or  by  doing  the  maximum  number  of 
iterations,  FNFIT  restarts  by  looking  for  a new  control  card. 

DGUK(X,  J,  DG)  - Derivati'^es 

DGDK  takes  the  derivative  of  the  function  F(X;  A(l),  A(2)...)  with 
respect  to  A(J)  at  X by  fitting  a 5th  order  polynomial  to  F evaluated 
at  5 values  of  A(J) . For  J = 0,  the  derivative  is  taken  with  respect 
to  X.  The  five  A(J)  values  are  A(J),  A(J)  ^ lo"^  C(J),  and  A(J)  + lO"^ 
C(J),  where  the  C(J)  are  the  "magnitudes"  of  the  A(J)  which  are  user 
supplied  (refer  to  USERS  MANUAL  section).  The  big  spread  in  the  points 
is  an  attempt  to  fit  both  functions  that  change  rapidly  and 
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functions  that  change  slowly. 

Having  fit  F(X;  A(l),  A(2)...)  by  a^  A(J)^  + a^  ACJ)*^  + A(J)^  + 

a^  A(J)  + a^,  the  derivative  is  given  by 

DF/DA(J)  = 4a^A(J)^  + 3ajA(J)^  + 2a2A(J)  + a^ 

GRAD  - Steepest  Gradient  Search 

GRAD,  like  FNFIT,  sets  up  the  array,  KL(J) , of  variable  parameter 
indicies.  Hence,  A(J)  refers  to  NV  actually  varied  parameters. 

2 

The  idea  of  GRAD  is  to  detect  the  steepest  gradient  in  the  X 
surface,  take  a step  down  that  gradient,  and  recompute.  The  size  of 
the  step,  STPSZ,  is  chosen  internally  as  0.4  in  the  non-dimensionalized 
parameter  space  described  below. 

A difficulty  in  this  kind  of  search  is  treating  all  parameters 

equally.  A step  size  of  10m  is  huge  for  a parameter  whose  range  is 

from  -lO’  to  10"^m,  is  negligible  for  one  whose  value  is  in  the 

order  of  lO^m,  and  is  meaningless  to  one  whose  dimension  is  grams. 

2 

Therefore,  GRAD  converts  the  problem  from  A-space  (X  , A(l),  A(2)...) 
to  B-space  (x^,  B(l) . B(2)...),  where  B(J)  = A(J)/C(J).  C(J)  refers, 
as  in  DGDK,  to  the  user  supplio'l  "magnitudes*'  of  the  A(J).  If  the 
C(J)s  fairly  well  describe  the  magnitudes  of  the  A(J)s,  then  a step 
of  1 in  B-spacc  has  the  same  relative  effect  on  all  the  parameters. 

Having  converted  to  B-space,  GRAD  proceeds  as  follows; 

2 

1)  like  FNFIT,  it  evaluates  G and  X and  outputs  as  directed.  GRAD 

exits  here  a)  if  MXSTP  steps  liave  been  taken  or  b)  if  two  successive 

2 

steps  have  resulted  in  a worse  X . Upon  any  exit,  the  values  of  the 
A(J)  are  set  to  those  that  gave  the  smallest  X^. 

1 

2)  the  next  step  is  *o  evaluate  thedX“/dB(Jl . 

Since 

<?X^  _ dX^  . «iA(J)  = dX^ 

aB(J)  " SA(JT  dB(J)  aA(J) 


• C(J) 


LvMfcy — "■ — ...  i^  ,.  ..^..jiiwtp  ^ _ 1 {.."•u-_  ;i»K4»i'i' -!?» I'S" 


2 2 
GRAD  can  evaluate  3x  /3B(J)  by  calling  UGDK  and  evaluating  3x  /3A(J)  as 

does  FNFIT.  Also,  the  root-square,  ( (3x‘’/3B(d))  ) ^ > is  computed. 

3)  The  step  in  B-space  is  taken  by  adding  to  each  B(J)  a change 

2 

equal  to  -(3x  /3B(.J))  • (STPSZ) /root-square.  Since  the  B(J)j  are  orthogonal, 
it  is  seen  that  this  procedure  takes  a step  in  B-space  of  size  STPSZ,  with 
the  largest  components  of  the  step  in  the  B(J)s  having  the  steepest 
gradient.  The  minus  sign  assures  movement  down  the  gradient. 

4)  GRAD  now  goes  about  checking  the  limits  and  second  derivative 

of  the  new  B(J)s.  First  COOLIT  is  called  to  check  maxima  and  minima  (since 
COOLIT  works  in  A-space,  the  A array  is  first  updated).  Then,  the  most 
recent  3x*'/3B(J)s  are  compared  to  the  previous  3x  /3B(J)s.  If  any 
parameter  has  either  exceeded  a COOLIT  limit  or  found  a shallower 
gradient,  that  parameter  is  hold  fixed,  and  the  program  returns  to 
step  3)  and  recalculates  the  B-step.  This  procedure  is  continued 
until  a)  an  allowed  step  is  taken,  or  until  b)  all  parameters  have  found 
a shallower  gradient  and  GRAD  exits.  If  a step  is  taken,  GRAD  loops 
back  to  step  1. 

COOLIT  - Parameter  Limiting  Subroutine 

As  described  in  the  USERS  MANUAL,  the  user  can  input  limits  on  the 
minimum  or  maximum  that  a parameter  can  attain,  or  limit  the  amount  a 
parameter  can  cliange  in  any  FNFIT  iteration.  (This  latter  limit  keeps 
the  search  from  making  wild  jumps  because  of  linear  approximations  to 
touchy  parameters.)  IVlien  COOLIT  is  called,  it  checks  the  A(J)s,  or  DA(o)s, 
versus  their  limits.  If  any  limits  are  exceeded,  the  A(J)  (or  DA(.J))  is 
reduced  to  its  allowed  limit.  Then,  the  errant  A(J)  is  placed  at  the 
top  of  the  list  of  parameters  temporarily  held  constant  (not  changed  in 
the  succeeding  iterations  in  FNFIT  or  GRAD). 

The  list  is  maintained  in  FNFIT  by  these  rules. 

1)  If  any  parameter  was  added  to  the  list  during  the  past  iteration 
no  parameter  is  released. 

2)  If  no  parameter  was  added  to  the  list  during  the  past  iteration 
then  the  parameter  at  the  bottom  of  the  list  (most  time  on  the  list)  is 


released. 


i- 


I 


r 

h 


s 

* 

E 


The  result  of  this  procedure  is  to  allow  the  program  to  readjust 
to  any  "artificial"  changes  caused  by  imposing  limits. 

To  avoid  the  possibility  of  a subtle  kind  of  loop,  the  limit  on 
the  DA(J)s  is  reduced  as  a function  of  the  iteration.  This  is  not  a 
forced  convergence;  it  merely  assures  non-repetition. 

IV.  GENERAL  CAUTIONS 

SQ  is  a function  defined  on  an  NPARAM  dimensional  space.  Depending 
on  the  form  of  G and  on  the  data,  there  may  exist  any  number  of  "local" 
minima,  i.e.,  sets  of  values  for  the  parameters  for  which  any  small 
change  in  any  of  the  parameters  produces  an  increase  in  SQ.  It  is 
assumed  that  the  best  fit,  and  hence  the  desired  set  of  parameters, 
corresponds  to  only  one  of  these  local  minima.  The  local  minimum  which 
the  program  seeks  out  depends  on  the  starting  set  of  parameters.  There 
exists  no  completely  satisfactory  method  to  avoid  the  multiple  minima 
problem;  it  is,  in  fact,  another  aspect  of  the  generally  unresolved 
problem  of  unambiguously  defining  a best  fit.  Two  features  of  FNFIT 
aid  in  best  fit  selection:  the  output  subroutine  displays  data  versus 

F(X)  for  easy  comparison  by  the  user,  and  initial  parameters  are  easily 
changed  to  "map  out"  local  minima. 

.Another  plienomenon  of  this  kind  of  fitting  program  occurs  when  F(X) 
is  a much  stronger  function  of  some  parameters  than  others.  In  such  a 
case,  some  parameters  are  varied  much  less  than  others.  The  hold- 
parameters-fixed  and  check-new-values-bei ore-re iteration  features  of 
FNFIT  are  useful  aids  in  assuring  a suitable  search  on  all  parameters. 

An  error  to  bo  avoided  in  any  function  fitting  problem  is  an 
attempt  to  use  indeterminant  parameters.  For  example,  the  function 
A(1)*(A(2)*X+A(3))  can  never  be  uniquely  fit,  since  A(2)  and  A{5) 
can  compensate  for  any  value  of  A(l).  Such  errors  can  be  very  subtle 
(e.g.  a variable  starting  point  for  a Fourier  Series  is  indeterminant.) 

A bad  start  can  also  fool  the  program  into  thinking  that  such  an  error 
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exists.  For  example,  in  fitting  a function  like  Aj-A2C  '3  (See  section 
V) , a search  that  starts  by  making  A^  and  A2  very  large  might  result  in 
values  for  Aj  that  are  so  small  that  the  function  is  approximately 
Ai  - A2,  and  acts  indeterminant.  Therefore,  the  size  of  changes  (Ll=l 
option  on  page  10)  should  be  held  to  the  minimum  required. 

If  the  function  does  not  closely  fit  the  data,  a possible  symptom 

of  indeterminancy  in  FNFIT  is  non-conve?gence.  If  the  data  does  fit 

well,  a probable  sympton  is  the  detection  of  a singular  matrix  by 

Subroutine  LSIMEQ.  If,  during  the  matrix  reduction- inversion  process, 

-12 

the  largest  pivot  eleraent  is  smaller  than  EPS1(10  , as  listed)  LSIMEQ 

aborts  and  returns  to  FNFIT.  FNFIT  prints  out  the  partially  reduced 
matrix,  the  indicies  of  the  previous  pivot  elements,  and  then  looks  for 
new  data  to  begin  a new  fit. 


V.  A FITTING  EXAMPLE 

As  an  example  of  the  operation  of  the  program,  twenty-five  values 
of  X(I),  (0  ^ X ^ 10.),  and  the  corresponding  values  of  Y(I)  from  the 
relationship  Y(I)  =1.  - EXP  (-X(I))  were  fit  with  the  following 

functions. 

A.  G = A(l)  ♦ A(2)  * EXP  C-A(3)  * X)  A(4)  * EXP  (A(5)  * X) 

4 

B.  G A(J)  * SIN  (J  * A(5)  * X) 

J=] 

5 

C.  G A(J)  * P (X) 

J=1 

where  P^fX)  is  the  Legendre  polynomial 

D.  G = A(l)  + A(2)  * ALOG  ( A(3)  * X + 1)  ♦ A(4)  * ALOG  ( A(5)  * (X**2)  + 1) 

Two  attempts,  with  different  starting  values  and  limits,  were  made 
with  case  B. 

The  sets  of  cards,  with  running  explanation  of  the  function  of  each, 
are  included  in  Appendix  B. 

The  same  set  of  starting  parameters  was  used  in  all  cases 


(A(l)  thru  A(4j  = 1.0,  A(5)  = 0.02),  indicating  the  relative  insensitivity 
of  the  search  in  this  case  to  initial  values.  Each  starting  value 
also  served  as  its  magnitude  ($  option  in  Set  B) . The  maximum  change 
in  any  iteration  was  limited  to  1/2  magnitude  ($  option  in  Set  C) . 

Also,  parameters  A(3)  and  A(5)  were  forced  to  stay  positive  in  cases 
A and  D,  and  also  in  case  B for  A(5) . 

The  initial  searches  were  limited  to  200  iterations.  However, 
cases  A and  B had  not  finished  (i.e.  exit  was  by  maximum  iterations 
rather  than  by  non-changing  values.)  Case  A is  difficult  in  that  A(4)  = 
0.0  is  the  exact  solution;  however,  as  A(4)  approaches  0.0,  A(5)  becomes 
indeterrainant.  Case  B was  then  continued  for  600  iterations  and  did  not 
improve;  the  search  was  in  a very  slowly  damped  oscillation  about  the 
final  values. 

The  results  are  plotted  in  Figure  2. 

The  values  of  for  the  cases  were: 


2.139  X 10 


6.903  X 10 


6.692  X 10 


4.951  X 10 
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Figure  2 * Results  of  Fitting  Example 
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APPENDIX  A 


I HEADING 

I SUBROUTINE  HEADING(  IK  ) 

I COHMON/PARAfl/NPARAM,  A(20),  SQ 

I INTEGER  V(7) 

I IF(  IK  .GT.  0 ) GO  TO  40 

s ■ NNN  « NPARAM 

I IF(  NNN  .GT.  10  ) NNN=10 

I ENCODE  ( S79,  V ) NNN 

I • 579  FORMAT  ( '(SHOIT*  . 'ER,5X,'  , *2HSQ,  13,  »(6X2HA'  , '(,12,1'  , 

I $ 'H)  ) )•  ) 

I 40  WRITE  ( 6,  V ) ( I,  I = 1,  NNN  ) 

I RETURN 

I END 


I FN 

i 

I SUBROUTINE  FN 

I COMMON/FUNCT/  F(IOOO) 

COMMON/DATA/NDATA,X(1000) ,Y(1000) ,W(1000) 
I DOUBLE  PRECISION  G,  T 

DO  70  1=1,  NDATA 
T + X(I) 

! CALL  FCN(  T,  0,  0.,  G ) 

i 70  F(I)  = G 

RETURN 

I END 

£ 

t 

I 

¥ 

I 

i DFN 

S 

I 

I SUBROUTINE  DFN (2) 

I COMMON/PARAM/NPARAM,  A(20),  SQ 

I COMMON/  DFUNCT/  DF(20),  NV,  K» (20) 

I DOUBLE  PRECISION  T 

I T = Z 

I DO  100  K = I,  NV 

5 ■ J = KLlK) 

! CALL  DGDK(  T,  J,  DG  ) 

J 100  DF(J)  = DG 

‘ RETURN 

^ END 


COOLIT 
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SUBROUTINE  COOLIT(  ITER,  KEY,  MARK  ) 

COMMON/  HOLD  / NHOLD,  IH0LD(20) 

COMMON/PARAM/NPARAM,  A(20),  SQ 

COMMON/  COOL/  NT,  MAXIT,  DA(20),  NDAM(3),  JDAM(3,20),  DAM(3,20) 
IF(  MARK  ) 100 

C CHECK  AMINS 
N2  = NDAM(2) 

IF(N2)  40,40 
DO  30  J2  = 1,  N2 
K2  = JDAM(2,  J2) 

IFf  A(K2)  .GE.  DAM(  2,  J2  ) ) GO  TO  30 
A(K2)  = DAM(2,  J2) 

KEY  = 1 
IF (NT)  28,20 
N1  =:>  NHOLD  + 1 
NN  = NHOLD  + NT 
DO  27  KK  =NI,NN 

27  IF(  IHOLD(KK)  .EQ.  K2  ) GO  TO  30 

28  NT  = NT+1 

IHOLD(  NHOLD  + NT  ) = K2 
30  CONTINUE 
C CHECK  AMAXS 

40  N3  = NDAM(3) 

IF(  N3  .EQ.  0 ) RETURN 
DO  60  J3  = 1,  N3 
K3  = JDAM(3,  J3) 

IF(  A(K3)  .LE.  DAM(3,  J3)  ) GO  TO  60 
A(K3)  = DAM(3,  J3) 

KEY  = 1 
IF(NT)  48,  48 
N1  = NHOLD  + 1 
NN  = NHOI.D  + NT 
DO  47  KK  =N1,  NN 

47  IF(  IHOLD(KK)  .EQ.  K3  ) GO  TO  60 

48  NT  = NT+1 

IHOLD(  NHOLD  + NT  ) = K3 
60  CONTINUE 
RETURN 

C CHECK  DA  S 
100  N1  = NDAM(l) 

IF(  N1  .EQ.  0 ) RETURN 

PON  = FLOATC  iter  ) / FLOAT(  MAXIT  ) 

DO  130  J1  = 1,  N1 
K1  = JDAM(  1,  J1  ) 

YUK  * DAM(1,  Jl)  + EXP(  - PON  ) 

IF(  ABS(  DA(Kl)  ) .LT.  YUK  ) GO  10  130 
DA(Kl)  = SIGNC  YUK,  DA(Kl)  ) 

KEY  = 1 

IF(NT)  128,128 
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N1  = NHOLD  + 1 
NN  = NHOLD  + NT 
DO  127  KK  =N1,  NN 

127  IF(  IHOLD(KK)  .EQ.  K1  ) GO  TO  130 

128  NT  NT+1 

IHOLDC  NHOLD  + NT  ) = K1 
130  CONTINUE 
RETURN 
END 


GAB 


SUBROUTINE  GAB(  MXSTP,  MAXIT,  IPRT,  ICHK,  NHOLD,  K2K  ) 
DIMENSION  GAB ( 140  ) 

NEW  = 1 


10  READC  5,  9,  ERR=25,  END«30  ) MXSTP,  MAXIT,  IPRT,  ICHK,  NHOLD 
9 FORMAT ( ) 

IF(  K2K  ) 5,  5 

PRINT  89,  ( GAB(M),  M = 1,  K2K  ) 

89  FORMATC  •!•,  ( T20,  13A6,A2,  / ) ) 

PRINT  98 

99  FORMATC  13A6,  A2  ) 

98  FORMAT  ( //  ) 

S RETURN 
25  CONTINUE 

IF(  NEW  ) 35,  35 
K2K  = 14 

KIK  = 1 

GO  TO  40 

35  KIK  = KIK  + 14 

K2K  = K2K  + 14 

40  READC  0,  99  ) ( GABCM),  M = KIK,  K2K  ) 

NEW  = 0 

GO  TO  10 

30  STOP 

END 
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i DO-FNFIT 
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LIB  MISD*LIBRARY. 

IN  AMUCK*FNFIT.,TPF$.FCN 

END 
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FNFIT 

C FN.FIT  MAIN  PROGRAM 

COMMON/  COOL/  NT,  MAXIT,  DAC20),  NDAM(3),  JDAM(3,20),  DAM(3,20) 
COMMON/DATA/NDATA,X(1000) ,Y(1000) ,W(1000) 

COMMON/  DFUNCT/  DF(20),  NV,  KL(20) 

COMMON/FUNCT/  F(IOOO) 

COMMON/  GRAUC  / MXSTP,  IPRT,  NNN,  ICHK 
COMMON/  HOLD  / NHOLD,  IH0LD(20) 

COMMON/  MAG  / C(20) 

COMMON/PARAM/NPARAM,A(20) ,SQ 
DATA  EPSl,  EPS2  / l.OE-12,  l.OE-5  / 

LOGICAL  PRE,  POST 

DIMENSION  AM(22, 22),  IR(22) , JC(22) ,RDA(22) 

DIMENSION  MAG (20) 

DIMENSION  ISDF0(20) 

K2K  = 0 

3 CALL  GAB(  MXSTP,  MAXIT,  IPRT,  ICHK,  NHOLD,  K2K  ) 

IF(  NHOLD)  4,4 

READ  (5,9)  ( IHOLD(K),  K = 1,  NHOLD  ) 

4 CONTINUE 

READ(  5,  9,  ERR  = 7777,  END  = 211  ) NNNN 

NPARAM  = NNNN 

DO  17  J = 1,  NPARAM 

READ(  5,  9,  ERR  = 18  ) A(J),  MAG(J) 

C(J)  » 10,**MAG(J) 

GO  TO  17 

18  C(J)  = A!’,S(  A(.J)  ) 

17  CONTINUE 
21  CONTINUE 

DO  11  I :=  1,3 
11  NDAM(I)  = 0 
PRE  = .FALSE. 

POST  = .FALSE. 

IF(  ICHK  ) 5,5 
DO  13  J = 1,  60 

READ(  5,  9,  ERR=8888,  END=15  ) Jl,  LI,  DUM 
NDAM(Ll)  = NDAM(Ll)  + 1 
N1  = NDAM(Ll) 

JDAM(L1,  Nl)  = Jl 
DAW(L1,  Nl)  = DUM 
GO  TO  13 

IF  COOLIT  CARD  IS  $,*,  ETC.,  ALL  PARAMETER  CHANGES  ARE  LIMITED 
TO  0.5*MAGNITUDE  OF  THE  PARAMETER 
8888  DO  14  J2=  1,  NPARAM 
NDAM(l)  = NDAM(l)  + 1 
Nl  = NDAM(l) 

JDAM(  1,  Nl  ) = J2 
DAM(  1,  Nl  ) - .5*C(J2) 
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14  CONTINUE 

13  CONTINUE 

15  IF(  NDAM(l)  .GT.  0 ) PRE  = .TRUE. 

IF(  ( NDAM(2)  .GT.  0 J .OR.  ( NDAM(3)  .GT.  0 ) ) POST  = 
5 CONTINUE 

9 FORMAT  ( ) 

CALL  DATAIN 
C 

C SET  UP  DGDK 

CALL  DGDK(  0.0,  -1,  0.0  ) 
r 

C PRINT  HEADING  INFO 

IF  ( NHOLD  .EQ.  0 ) GO  TO  48 
PRINT  249,  ( IHOLD(L),  L = 1,  NHOLD) 

48  CONTINUE 

PRINT  609 

IF(  NDAM(l)  ) 620,  620 

NDAMl  = NDAM(l) 

PRINT  619,  ( JDAM(1,K),  DAM(1,K),  K =■•  1,  NDAMl  ) 

620  IFC  NDAM(2)  ) 630,  630 

NDAMl  = NDAM(2) 

PRINT  629,  ( JDAM(2,K),  DAM(2,K),  K = 1,  NDAMl  ) 

630  IF(  NDAM(3)  ) 640,  640 

NDAMl  = NDAM(3) 

PRINT  639,  ( JDAM(3,K),  DAM(3,K),  K = \ NDAMl  ) 

640  CONTINUE 

249  FORMAT  ( 29H0PARAMETERS  HELD  CONSTANT...  , 20 ( 14,  IH, 
609  FORMAT  ( ‘0»  ) 

619  FORMAT  ( ' CHANGE  LIMITED  FOR  PARAMETERS  ...  • 

$ 4(  T37,  5(  14,  » (',  1PE9.2,  •),')/)) 

629  FORMAT  ( ' MINIMUM  LIMITED  FOR  PAlL'vMETERS  ...  • 

$ 4(  T37,  5(  14,  ’ (',  1PE9.2,  '),')/)) 

639  FORMAT  ( ' MAXIMUM  LIMITED  FOR  PARAMETERS  ...  « 

$ 4(  T37,  5(  14,  ' (’,  1PE9.2,  ’),*)/)) 

C 

C CALL  GRADIENT  SEARCH 

NNN  = NDATA  + NHOLD  - NPARAM 
IF(  MXSTP  ) 40,  40 
CALL  GRAD 
40  PRINT  299 

299  FORMAT  ( / / ' CHISQ  MINIMIZATION  SEARCH  » ) 

CALL  HEADING (MXSTP) 

KEY  = 0 
NT  = 0 
C 

C THE  ITERATIVE  LOOP  STARTS  HERE 

C 

ITER  = 0 
N1  = NHOLD  + 1 
1002  CONTINUE 


.TRUE. 


) ) 
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NOW  CALCULATE  SQ 
SQ  = 0.0 

DO  20  I - 1,  NDATA 

10  SQ  = SQ  + W(I)*(  Y(I)  - F(I)  )**2  / NNN 

OUTPUT  VARIABLES  FROM  PREVIOUS  ITERATION 
CALL  OUTPUT ( ITER,  IPRT  ) 

1001  CONTINUE 

IF  NO  NEW  PARAMETER  IS  BEING  HELD  THIS  ITER.,  REMOVE  THE  OLDEST 
IF(  ( KEY  .NE.  0 ) .OR.  ( NT  .EQ.  0 ) ) GO  TO  140 
NT  = NT  - 1 

NAll  = NHOLD  + NT 
DO  135  KK  = Nl,  NAH 
135  IHOLD(KK)  = IH0LD(  KK+1  ) 

140  CONTINUE 

NAH  = NHOLD  + NT 
NV  = NPARAM  - NAH 

MAKE  SURE  THAT  SOME  PARAMETER  IS  ALLOWED  TO  VARY 
IF(  NV  ) 195,195,200 
195  ITER  = ITER  + 1 
PRINT  199,  ITER 

199  FORMATC  IX,  13,  * ALL  VARIABLES  HELD  • ) 

KEY  = 0 

GO  TO  1001 

200  CONTINUE 

IF  ( ITER  .EQ.  MAXIT  ) GO  TO  1010 

NOW  SET  UP  ARRAY,  KL,  OF  VARIABLE  PARAMETER  INDICES,  IN  NUMERICAL 
ORDER 
K = 1 

DO  210  II  = 1, NPARAM 

IF(  NAH  .EQ.  0 ) GO  TO  215 

DO  220  JJ  = 1,NAH 

220  IF(  IHOLD(JJ)  .EQ.  II  ) GO  TO  210 
215  KL(K)  = II 
K = K+l 
210  CONTINUE 

SET  UP  NV  X NV+1  MATRIX  OF  COEF.  , ONE  DATA  POINT  AT  A TIME 
NVl  = NV+1 
DO  50  II  = 1,NV 
ISDFO(Il)  = 0 
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DO  50  12  = 1,  NVl 
50  AM(I1,I2)  = 0. 

DO  100  I = l.NDATA 
CALL  DFN(X(I)) 

C INTO  / DEFUNCT/  APPEARETH  ( DF(  KL(J)  ) , J = 1,NV  ) EVALUATED  AT  X(I) 
D(-  80  K1  = 1,NV 
K = KL(Kl) 

C 

C CHECK  IF  DFCSOME  VARIABLE)  IS  I.D.  = 0 

IF(  ABS(  DF(K)  ) .LT.  EPSl  ) ISDFO(Kl)  = ISDFO(Kl)  + 1 
C 

C FILL  UP  COEFICIENT  MATTRIX  FOR  DA 

AM(K1,  NVl  ) = AMCKl,  NVl  ) ♦ W(I)*(  Y(I)  - F(I)  )*DF(K) 

DO  60  K2  = 1,NV 
L = KL(K2) 

AM(K1,K2)  = AM(K1,K2)  + W(I)*DF(K)*DF(L) 

60  CONTINUE 
80  CONTINUE 
100  CONTINUE 
C 

C FINISH  CHECK  FOR  DF  I.D.  = 0 AND  REDUCE  MATRICES  OF  VARIABLES 
C ACCORDINGLY 

NVQ  = NV 

DO  300  K = 1,  NVQ 

IF(  ISDFO(K)  .NE.  NDATA  ) GO  TO  300 

PRINT  599,  KL(K),  KL(K) 

599  FORMATC  • DF/DA(',  12,  •)  IS  I.D.  = 0 NEXT  ITER.  A(»,  12.  •)  WI 
$lL  BE  HELD.  ' ) 

NVl  = NV 
NV  * NV-1 
IF(NV)  195,  195 
DO  340  LI  = K,  NV 
KL(Ll)  = KL(  Ll+1  ) 

DO  340  L2  = 1,  NVl 
340  AM(  LI,  L2  ) = AM(  Ll+1,  L2  ) 

DO  350  L2  = K,  NVl 
DO  350  LI  = 1,  NV 
350  AiM(  LI,  L2  ) = AM(  LI,  L2+1  ) 

300  CONTINUE 
C 

C NOW  FIND  SOLN  TO  DA(20) 

CALL  LSIMEQ  (AM,  22,  IR,  JC,  NV,  EPS1,RDA,  lERRl  ) 

IF(IERR1  .EQ.  -1  ) GO  TO  7734 
DO  102  J = 1,  NPARAM 
102  DA(J)  =0.0 

DO  108  K1  = 1,  NV 
K2  = KL(Kl) 

108  DA(K2)  = RDA(Kl) 

C 

C QUIT  IF  NO  PARAMETERS  ARE  CHANGING  OR  ALL  ARE  ZERO 

IF(  ( ITER  .LT.  4 ) .OR.  ( t’l  .GT.  0 ) ) GO  TO  106 
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IF  ANY  PARAMETER  CHANGE  IS  LIMITED,  CALL  COOLIT 
IF(  PRE  ) CALL  COOLIT(  ITER,  KEY,  -1  ) 


NOW  ADD  THE  DA  S TO  THE  PROPER  A S 
110  DO  120  J = 1,  NPARAM 
120  A(J)  = A(J)  + DA(J) 

IF  ANY  PARAMETER  IS  LIMITED,  CALL  COOLIT 
IF(  POST  ) CALL  COOLIT ( ITER,  KEY,  1 ) 

1000  CONTINUE 

ITER  = ITER  + 1 

IF(  ITER  .LE.  MAXIT  ) GO  TO  1002 

THE  ITERATIVE  LOOP  ENDS  HERE. 

1010  CALL  OUTPUTC  ITER, 4 ) 

GO  TO  3 

7734  PRINT  689, 

689  FORMAT  C * LARGEST  ELEMENT  .LT.  EPSl  IN  LSIMEQ  ' ) 

DO  701  I = 1,NV 

PRINT  9,  ( AM(I),J),  J = 1,  NVl  ) 

701  CONTINUE 
PRINT  797 

797  FORMAT  ( ' VARIABLES  THIS  ITERATION  ' ) 

PRINT  799,  ( KL(J),  J = 1,  NV  ) 

PRINT  798 

798  FORMAT(  ‘OINDETERMINATE****  PROBLEM  MIGHT  BE..,',/, 

$ T20,  'IJ  PARAMETERS  ARE  NOT  INDEPENDENT',/, 

$ T20,  '2)  NUMBER  PARAMETERS  + 1 .GE.  NUMBER  DATA  POINTS',  /, 

$ T20,  '5)  DATA  NOT  IN  RIGHT  RANGE  TO  DETERMINE  ONE  PARAMETER'  ) 
C 

799  FORMAT  ( / IX  20IS  ) 

GO  TO  1010 

211  DO  212  J = 1,  NPARAM 

212  A(J)  = C(J) 

GO  TO  21 

7777  PRINT  9999 

9999  FORMAT ( ' MESS-UP  IN  MAIN  PROGRAM  READ  ' ) 

STOP 

END 
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50 

52 


60 


DATAIN 

SUBROUTINE  DATAIN 

COMMON/DATA/NDATA,X(1000),y(1000),W(1000) 
C0»WN/  XTREME/  DELYIO,  DELXIO 
READ(  5,  9,  ERR=40,  END=»25  ) tYT 
AW  = ABS(  WT) 

DO  20  1=1,  1000 

READ(  5,  9,  ERR  = 7734,  END  = 33  ) XX,  YY 

X(I)  = XX 

Y(I)  = YY 

W(I)  = 1.0 

NDATA  = I 

CONTINUE 

CONTINUE 

IF(  AWT  .LT.  l.OE-8  ) GO  TO  52 
DO  50  1=1,  NDATA 

T = ABS(  Y(I)  ) 

IF(  T .LT.  l.OE-12  ) GO  TO  50 

W(I)  = T**W 

CONTINUE 

CONTINUE 

XMIN  = X(l) 

XMAX  = X(l) 

WORM  =0.0 

DO  60  1=1,  NDATA 

WORM  = WORM  + W(I) 

XMIN  = AMIN1(  XMIN,  X(I)  ) 

XMAX  = AMAX1(  XMAX,  X(I)  ) 

CONTINUE 

DELXIO  = ( XMAX  - XMIN  ) / ( 10.  * NDATA  ) 
DO  30  1=1,  NDATA 

W(I)  = NDATA*W(I)/  WORM 
RETURN 


MESS-UP  IN  DATAIN  READ  ' ) 


30 
25 

7734  PRINT  69 
STOP 

FORMAT  ( ) 

FORMAT ( 

CONTINUE 

DO  45  1 = 1,  1000 

READC  5,  9,  ERR  = 7734,  END  = 52  ) 
NDATA  = I 
CONTINUE 
GO  TO  7734 
END 


9 

69 

40 


45 


X(I)  , Y(I)  , W(I) 
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DGDK 


30 

32 


36 

34 

40 


10 

20 


SUBROUTINE  DGDK(  T,  K5,  DG  ) 
C0MMON/PARAM/NPARAM,A(20) ,SQ 
COJWON/  XTREME/  DELYIO,  DELXIO 
C0W10N/  MAG  / C(20) 

DIMENSION  SGN(2)/  -1.,  1.  / 

DOUBLE  PRECISION  G,  T,  DEL 

DOUBLE  PRECISION  P(2),  Q(3),  DT(21,3), 

DOUBLE  PRECISIOI  D3(21,3  ) 

IF(K5)  30,  40,  40 
DT(1,1)  = DELXIO 
DO  32  J = 1,  NPARAM 
DT(  J+1,  1 ) = 0.01*C(J) 

NPl  = NPARAM  + 1 
DO  34  N = 1,  NPl 
DT(N,2)  = DT(N,l)*1.0E-4 
DO  36  II  = 1,2 
D2(N,II)  = DT(N,II)**2 
D3(N,II)  = DT(N,II)**3 
CCX4TINUE 


D2(21,3),  DN0M(21) 


DNOM(N)  = 0.5/  C DT(N,1)*DT(N,2)*(  D2(N,2)  - D2(N,1'  ) ) 

CONTINUE 

RETURN 

CONTINUE 


N = KS  ♦ 1 
DO  20  L = 1,2 
DO  10  K = 1,2 
DEL  = DT(N,L)*SGN(K) 

CALL  FCN(  T,  K5,  DEL,  G ) 

P(K)  = G 

Q(L)  = P(2)  - P(l) 

DG  = ( Q(1)*D3(N,2)  - Q(2)*D3(N,1)  ) * DNOM(N) 

RETURN 

END 


OUTPUT 

SUBROUTINE  OUTPUT ( ITER,KKEY  ) 
C0MM0N/PARAM/NPARAM,A(20) ,SQ 
C0t«0N/FUNCT/F(1000) 

CO^WON/DATA/NDATA,X(1000) ,Y(1000) ,W(1000) 
DOUBLE  PRECISION  T 
DIMENSION  DD(IOOO),  BEST(20) 

LOGICAL  YEP 
YEP  = .TRUE. 

IF(  KKEY  ) 15,15 

32 
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695 


20 


25 

18 


35 

499 

15 

16 


IFC  KKEY  .NE.  4 ) GO  TO  18 
PRINT  695 

FORMATC  ///,  lOX,  'BEST  FIT’  / ) 

SQ  = SOFAR 
ITER  = I BEST 
DO  2U  J = 1,  NPARA.M 
A(J)  = BEST(J) 

CONTINUE 
CALL  FN 

DO  25  1=1,  NDATA 

T = X(I) 

CALL  DGDK(  T,  0,  D ) 

DD(I)  = D 

CONTINUE 

CONTINUE 

IF(  NPARAM  .GT.  10  ) GO  TO  30 

PRINT  499,  ITER,  SQ,  ( A(I),  I = 1,  NPARAM  ) 

CONTINUE 

FORMAT  ( IX,  13,  IX,  11(1PE11.3)  ) 

IF(  ( ITER  .EQ.  0 ) .OR.  C SQ  .LT.  SOFAR  ) ) 
CONTINUE 

IF(  KKEY  .LE.  1 ) RETURN 
IF(  K^’EY  .GE.  3 ) GO  TO  6 
IF(  YEP  ) RETURN 


GO  TO  12 


6 PRINT  598 

598  FORMATC  /,54X, ' I ' ,5X, ’X(I) ' , 8X,  ’Y(I)’,  8X,  ’FCD’,  8X,  ’W(I)’ 
IF (KKEY  .NE.  4 ) GO  TO  8 

PRINT  698 

698  FORMAT  ( ’+'  , T107,  ’DF/DX(I) ’ ) 

8 PRINT  697 

697  FORMATC  //  ) 

IF (KKEY  .EQ.  4 ) CO  TO  9 

PRINT  599,  C I,  XCI),  Y(I),  F(I),  K(I),  1=1,  NDATA  ) 

599  FORMAT  ( / , 1000(50X,I4,4C2X1PE10.4)  / ) ) 


GO  TO  7 

9 CONTINUE  , , , 

PRINT  699,  C I,  X(I),  Y(I),  FCD,  W(I),  DDCI),  I = 1, NDATA  ) 
699  FORMAT  ( / , 1000(50X, I4.5C2X1PE10.4)  / ) ) 

7 CONTINUE 

RETURN 

12  SOFAR  = SQ 
I BEST  = ITER 


DO  14  J = 1,  NPARAT-I 
BEST(J)  = A(J) 

14  CONTINUE 

YEP  = .FALSE. 


GO  TO  16 
30  CONTINUE 

PRINT  499,  ITER,  SQ,  CA(I),  I = 1,10  ) 
PRINT  799,  (ACI),  I = H,  NPARAM  ) 

799  FORMATC  T15,  10C1PE11.3)  ) 

GO  TO  35 
END 


) 
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SUBROUTINE  FCN(  X,  KS,TDA,  G ) 

COMMON/PARAM/NPARAM,  A(20),  SQ 

DOUBLE  PRECISION  G,  X , TDA 

IF(  K5  ,EQ.  0 ) X = X +TDA 

IF(  K5  .GT.  0 ) A(K5)  = A(K5)  +TDA 

CONTINUE 

INSERT  FUNCTION  ( G = FUNCTION(  X,  A(1),A(2),...  ) HERE  BETWEEN 
TOE  CONTINUES 
CONTINUE 

IF(  K5  .EQ.  0 ) X = X - TDA 

IF(  K5  .GT.  0 ) A(K5)  = A(KS)  - TDA 

RETURN 

END 


GRAD 

SUBROUTINE  GRAD 

COKWON/DATA/NDATA,X(1000) ,Y(1000) ,W(1000) 

CO^WON/  COOL/  NT,  MAXIT,  DA(20),  NDAM(3),  JDAM(3,20),  DAM(3,20) 
COMMON/  DFUNCT/  DF(20),  NV,  KL(20) 

COMMON/ FUNCT/  F(IOOO) 

CONWON/  GRADC  / MXSTP,  IPRT,  NNN,  ICHK 
COMMON/  HOLD  / NHOLD,  IHOLD(20) 

COMMON/ PAkAM/NPARAM, A( 20) , SQ 
COMMON/  MAG  / C(20) 

DATA  STPSZ  / 0.4  / 

DIMENSION  B(20),  D(20),  DXDB(20),  DXDB1(20),  ISG(20) 

LOGICAL  LI,  L2,  L3 

SET  UP  ARRAY,  KL,  OF  VARIABLE  PARAMETER  INDICES 
NV  = NPARAM  - NHOLD 
K = 1 

DO  100  J1  = 1,  NPARAM 
DXDB1(J])  = 0.0 
ISG(Jl)  = 1 
IF(  NHOLD  ) 90,  90 
DO  80  J2  = 1,  NHOLD 
80  IF(  IHOLD(J2)  .EQ.  J1  ) GO  TO  100 
90  KL(K)  = J1 
K = K + 1 
100  CONTINUE 
110  CONTINUE 

NHOLDl  = NHOLD  ♦ 1 
C 

C SET  UP  ARRAY,  B(J),  OF  DIMENSIONLESS  PARAMETERS 
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DO  150  J = 1,  NPARAM 
B(J)  = A(J)  / C(J) 

150  CONTINUE 
PRINT  299 

299  FORMAT  ( //  • STEEPEST  DESCENT  SEARCH  » / ) 

CALL  HEADNGC  0 ) 

>WM  = MXSTP  + 1 
LI  = .TRUE. 

L3  = .FALSE. 

C 

C THE  ITERATIVE  LOOP  STARTS  HERE 
C 

DO  1000  ISTP  = 1,  ^M( 

ITER  = ISTP  - 1 
NT  = 0 
NF  = NHOLD 

DO  41  J = NH0LD1«  NPARAM 
41  IHOLD(J)  = 0 
CALL  FN 
SQ  = 0. 

DO  20  1=1,  NDATA 

20  SQ  = SQ  + W(I)*(  Y(I)  - F(I)  )**2  / NNN 

C 

CALL  0UTPUT(  ITER,  IPRT  ) 

C 

IF(  SQ  .LT.  SOFAR  ) LI  = .TRUE. 

IF(  .NOT.  LI  ) GO  TO  56 
SOFAR  = SQ 

DO  50  J = I,  NPARAM 
50  D(J)  = A(J) 

GO  TO  560 
56  CONTINUE 

DO  550  J = 1,  NPARAM 
IF(  ISG(J)  .EQ.  -1)  GO  TO  555 
550  CONTINUE 
GO  TO  560 
555  PRINT  399 

399  FORMATC  » CHI-SQ  INCREASED  AND  GRADIENT  CHANGED  SIGN.  IF  NO  FINAL 
$ FIT,  CHECK  USER  SPECIFIED  MAGNITUDES.  » ) 

560  CONTINUE 

IF(L3)  GO  TO  21 

IF(  LI  .OR.  L2  ) GO  TO  25 

21  CONTINUE 

DO  35  J = 1,  NPARAM 
35  A(J)  = D(J) 

RETURN 

25  CONTINUE 

IF(  ISTP  .EQ.  ^•W  ) RETURN 
L2  = LI 
LI  = .FALSE. 
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30 


300 

400 

401 


500 


FIND  TOE  PARTIAL  DERIVATIVES  OF  CHISQ  W.R.T.  TOE  B(J)  S 
DO  30  J = 1,  NPARAM 
DXDBCJ)  = 0. 

DO  400  1=1.  NDATA 

DFN  FILLS  DFCi;; 

CALL  DFN(X(I)) 

DO  300  K1  = 1,  NV 
K2  = KL(Kl) 

DXDB(K2)  = DXDBCK2)  + W(I)*DF(K2)*C(K2)*(  F(I)  - Y(I)  ) 

CONTINUE 

CONTINUE 

DXB  = 0. 

DO  500  K1  = 1,  NV 
K2  = KL(Kl) 

DDD  = DXDB(K2)*DXDB1(K2) 

ISG(K2)  = ISIGNC  1,  DDD  ) 

IF(  ABS(  DDD  ) .LT.  l.OE-24  ) ISG(K2)  = 1 
DXB  = DXB  + DXDB(K2)**2 
DXB  = SQRT(DXB  ) 


BOO 
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NOW  GET  NEW  COEFICIENTS 
DO  800  K1  = 1,  NV 
K2  = KL(Kl) 

B(K2)  = B(K2)  - DXDBCK2) 
CONTINUE 

DO  40  J = 1,  NPARAM 
A(J)  = B(J)*C(J) 
CONTINUE 
KEY  = 0 


* STPSZ  / DXB 


HOLD  ANY  PARAMETER  TOAT  EXCEEDS  ITS  »COOLIT»  LIMITS,  OR  WHO  IS 
' IN  THE  VALLEY  ' ( I.E.  WHOSE  DXDB  HAS  DECREASED  ) 

IF(  ICHK  ) 501,  501 
CALL  COOLITC  ITER,  KEY,  1 ) 

501  NF  = NHOLD  + NT 

DO  505  K1  = 1,  NV 
K2  = KL(Kl) 

DO  502  Ml  = NHOLDl,  NF 

502  IF(  IHOLD(Ml)  .EQ.  K2  ) GO  TO  505 

IF(  ABS(  DXDBCK2)  ) ,GE.  ABS(  DXDB1(K2)  ) ) GO  TO  504 
NT  = NT  + 1 
KEY  = 1 

IHOLDC  NHOLD  + NT  ) = K2 
GO  TO  505 

504  CONTINUE 
DXDB1(K2)  = DXDB(K2) 

505  CONTINUE 

IF(  KEY  .EQ.  0 ) GO  TO  1000 
NF  = NHOLD  + NT 
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IF(  NPARAM  - NF  ) 60,  60 

DO  42  J = NHOLDl,  NF 
K = IHOLD(J) 

B(K)  = A(K)/C(K) 

42  DXDB(K)  = 0. 

GO  TO  401 
60  L3  = .TRUE. 

1000  CONTINUE 
RETURN 
END 
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APPENDIX  B 

The  following  cards  do  case  A,  and  case  B twice,  as  described 
in  the  main  text. 


0ASG,A  AMUCK  * FNFIT. 

0ADD  AMUCK  * FNFIT. DO-G 

G=A(1)  + A(2)*EXP(-A(3)*X)  + A(4)*EXP(A(5)*X) 
0ADD  AMUCK  * FNFIT. DO-FNFIT 
G=A(1)  + A(2)*EXP(-A(3)*X)  + A(4)*EXP(A(5) *X) 

100,  200,  1,  1,  0 

5 

1.0,  $ 

1.0,  $ 

1.0,  $ 

1.0,  $ 

2.0E-2,  1 

$ 

3,  2,  l.OE-3 
5,  2,  l.OE-3 
f!£0F 
0.0 

0.0,  0.0 
0.1,  9.516E-2 
0.2,  1.826E-1 
0.3,  2.591SE-1 
etc. 

0EOF 

0EOF 


"GAB"  Card 

Control 

NPARAM 

initial 
values  and 
"magnitudes" 

COOLIT  Option 
Limit  A(3)5A(S) 
Greater  than  10 
Ends  COOLIT  Data 
WT  (0.0  means 
unweighted) 

Data  (X(I),  Y(I)) 


Ends  Data 
Ends  Program 


NOTE:  Must  end  program  bec’use  new  case  requires  new  absolute  element 

to  be  compiled. 


0ADD  AMUCK  * FNFIT. DO-G 
G = 0.0 

DO  10  J = 1,  4 

10  G = G + A(J)  * SIN  (J*A(5)*X) 
0ADD  AMUCK  * FNFIT. DO-FNFIT 


PrecediRg  page  Uank 
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THIS  IS  A FOUR  SINE  FIT 
G = SUM  (A(J).SIN(A(5).J.X)) 
100,  200,  1,  1,  0 
5 

1.0,  $ 

1.0,  $ 

1.0,  $ 

1.0,  $ 

0.02,  $ 

$ 

5,  1,  0.005 
5,  2,  0.001 
@EOF 
0.0 

\ 

0.0,  0.0 
0.1,  9.516E-2 
etc. 

@EOF 

100,  400,  1,  1,  0 


A(5) 


•'GAB" 

Cards 

Control 

:;?ARAM 

initial 

values 

and 

magnitudes 
COOLIT  Option 
Small  Change 
A(5)  > 10"^ 

End  COOLIT  Data 
WT 

DATA 


Control  for 
second  run 
Instead  of  NPARAM 
Repeats  w/old  "MAGS" 
COOLIT  Data 


5,  1,  0.005 
5,  2,  0.001 
@E0F 

j Footnote 


instead  of  weight 


0.0,  0.0,  1.0 

0.1,  9.516E-2,  1.1 

1 read  in  weights 

0.2,  1.826E-1,  1.3 

etc. 

j X(I),  Y(I),  W(I) 

@E0F 

§E0F 

End  Program 
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GLOSSARY  OF  INPUT  TERMS 


A(l),  A(2),  ... 
DA(D... 

DUM 

G 

ICHK 

IHOLD(i),  IHOLUC2)... 

IPRT 

J1 

LI 

MAG 

MAXIT 

MXSTP 

NHOLD 

NPARAM 

W(l),  W(2)... 

WT 


Starting  values  for  the  parameters 
Change  in  a parameter  in  an  iteration 
Actual  Value  of  a COOLIT  limit 
User  supplied  function 
COOLIT  indicator 

Indicies  of  paraiiieters  held  constant 

OUTPUT  indicator 

Index  of  parameter  to  be  cooled 

Type  of  COOLIT  limit 

Magnitude  (power  of  10)  of  a parameter 

Maximum  number  of  FNFIT  iterations 

Maximum  number  of  GRAD  steps 

Number  of  parameters  held  constant 

Total  number  of  parameters 

Weights  for  data  points,  if  input  by  user 

Power  of  Y(I)  to  weight  data  points,  if  computed 
by  program 


X(I)\ 

Y(I)J 


DATA  to  be  fit 


Precediiii  m*  M 
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